Mise en place

Chargement des packages:

library(lme4) # pour la modélisation multiniveau
library(multilevelTools) # pour l'analyse des modèles
library(lmerTest) # pour l'évaluation des modèles
library(sf) # pour les objets spatiaux (cartes)
library(tidyverse) # pour la manipulation des données
library(readxl) # pour ouvrir les fichiers xlsx
library(ggrepel) # pour des étiquettes de graphiques propres

Chargement des données:

Jeu de données de niveau 1

Il s’agit des vignobles de Bourgogne. On a dans un tableau les caractéristiques de 2391 vignobles en termes de:

prix moyen du vin (millésimes 1998, 2002 et 2003)
qualité pédologique et météo (radiations solaires, pluie, etc.)
appellation d’origine contrôlée en 5 niveaux (Côteaux bourgignons < Bourgogne < Village < Premier cru < Grand cru.)
surface

Sources: Ay J.S., Hilal M., 2020, « Les déterminants naturels et politiques des AOC viticoles de Côte-d’Or », Cybergeo: European Journal of Geography, document 973, DOI : https://doi.org/10.4000/cybergeo.36443 et https://www.hachette-vins.com/guide-vins/actualite-vin/487/les-meilleurs-millesimes-des-vins-de-france/

vignobles <- st_read("Data/Vignobles_prix.shp") %>% 
  st_transform(crs = "EPSG:2154")

## Reading layer `Vignobles_prix' from data source 
##   `/Users/ccottineau/Documents/GitHub/QuantiLilleGeo/Cours_VinBourgogne/Data/Vignobles_prix.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2391 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 826631.6 ymin: 6645988 xmax: 851465.7 ymax: 6690810
## Projected CRS: RGF93 v1 / Lambert-93

summary(vignobles)

##     Concat             LIBCOM              NOM               NIVEAU         
##  Length:2391        Length:2391        Length:2391        Length:2391       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##     SURFACE          SCORE_b         Source             PrixMoy       
##  Min.   :0.0000   Min.   :18.40   Length:2391        Min.   :    5.0  
##  1st Qu.:0.1300   1st Qu.:67.90   Class :character   1st Qu.:   23.0  
##  Median :0.3100   Median :73.20   Mode  :character   Median :   36.0  
##  Mean   :0.4773   Mean   :71.89                      Mean   :  154.4  
##  3rd Qu.:0.6100   3rd Qu.:76.60                      3rd Qu.:   87.0  
##  Max.   :7.2000   Max.   :94.22                      Max.   :25960.0  
##                                                      NA's   :4        
##           geometry   
##  MULTIPOLYGON :2391  
##  epsg:2154    :   0  
##  +proj=lcc ...:   0  
##                      
##                      
##                      
##

vignobles %>% 
  dplyr::filter(LIBCOM == "FLAGEY-ECHEZEAUX" & NIVEAU == "Premier cru")

Jeu de données de niveau 2

Il s’agit des communes de Côte d’Or. On a dans un tableau les caractéristiques de 31 communes en termes de:

superficie
population
côte d’appellation
hierarchie administrative

Source: https://geo.data.gouv.fr/en/datasets/cac9f2c0de2d3a0209af2080854b6f6a7ee3d9f4

communes <- st_read("Data/communes.shp",
                     stringsAsFactors = T) %>% 
  st_transform(crs = "EPSG:2154")

## Reading layer `Communes' from data source 
##   `/Users/ccottineau/Documents/GitHub/QuantiLilleGeo/Cours_VinBourgogne/Data/Communes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 31 features and 20 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 826428 ymin: 6645899 xmax: 855115 ymax: 6691133
## Projected CRS: RGF93 v1 / Lambert-93

summary(communes)

##       gid          id_geofla       code_comm    insee_com 
##  Min.   :  563   Min.   :  563   010    : 1   21010  : 1  
##  1st Qu.: 6744   1st Qu.: 6744   037    : 1   21037  : 1  
##  Median :14409   Median :14409   054    : 1   21054  : 1  
##  Mean   :11410   Mean   :11410   110    : 1   21110  : 1  
##  3rd Qu.:14602   3rd Qu.:14602   133    : 1   21133  : 1  
##  Max.   :14749   Max.   :14749   150    : 1   21150  : 1  
##                                  (Other):25   (Other):25  
##                  nom_comm               statut     x_chf_lieu     y_chf_lieu   
##  ALOXE-CORTON        : 1   Chef-lieu canton: 3   Min.   :8291   Min.   :66472  
##  AUXEY-DURESSES      : 1   Commune simple  :27   1st Qu.:8349   1st Qu.:66570  
##  BEAUNE              : 1   Sous-préfecture : 1   Median :8431   Median :66660  
##  BROCHON             : 1                         Mean   :8419   Mean   :66676  
##  CHAMBOLLE-MUSIGNY   : 1                         3rd Qu.:8486   3rd Qu.:66774  
##  CHASSAGNE-MONTRACHET: 1                         Max.   :8516   Max.   :66897  
##  (Other)             :25                                                       
##    x_centroid     y_centroid       z_moyen        superficie     population   
##  Min.   :8288   Min.   :66477   Min.   :222.0   Min.   :  90   Min.   : 0.20  
##  1st Qu.:8350   1st Qu.:66576   1st Qu.:254.0   1st Qu.: 687   1st Qu.: 0.35  
##  Median :8450   Median :66661   Median :293.0   Median : 891   Median : 0.60  
##  Mean   :8418   Mean   :66677   Mean   :303.6   Mean   :1147   Mean   : 2.21  
##  3rd Qu.:8482   3rd Qu.:66774   3rd Qu.:349.0   3rd Qu.:1267   3rd Qu.: 1.30  
##  Max.   :8514   Max.   :66899   Max.   :464.0   Max.   :3619   Max.   :21.90  
##                                                                               
##    code_cant code_arr code_dept      nom_dept  code_reg     nom_region
##  05     :8   1:23     21:31     COTE-D'OR:31   26:31    BOURGOGNE:31  
##  24     :7   2: 8                                                     
##  15     :6                                                            
##  23     :5                                                            
##  06     :2                                                            
##  38     :1                                                            
##  (Other):2                                                            
##             Cote             geometry 
##  Côte D'or    : 3   POLYGON      :31  
##  Côte de Beane:16   epsg:2154    : 0  
##  Côte de Nuit :12   +proj=lcc ...: 0  
##                                       
##                                       
##                                       
##

Visualisation des données:

Avec ggplot, utilisation de geom_sf pour les objets spatiaux.

ggplot(data = vignobles) +
  geom_sf(aes(fill = log(PrixMoy)), lwd=0.01, colour="white")  +
  scale_fill_viridis_c(option = "viridis")+
  labs(title = "Niveaux de qualité physique des vignobles")

Pour des étiquettes réactives, on utilise geom_text_repel:

ggplot(data = communes) +
  geom_sf(aes(fill = Cote), lwd=0.01, colour="white") +
  geom_text_repel(data = communes[communes$population > 1,],
                        aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Appartenance communale aux côtes bourguignonnes")

Création de variables niveau 2 à partir d’aggrégation de données niveau 1

Info sur la distribution des prix, des surfaces, des qualités physiques des vignobles, etc.

Vignobles_par_communes <- as_tibble(vignobles) %>%
  group_by(LIBCOM) %>%
  summarise(aire_vignoble = sum(SURFACE),
            MoyPond_b = sum(SCORE_b * SURFACE) / sum(SURFACE),
            N_Parcelles = n(),
            Prix_moyen = mean(PrixMoy,na.rm=T),
            Prix_median = median(PrixMoy,na.rm=T),
            )

Info sur les appellations d’origine contrôlée:

Vignobles_par_communes_par_AOC <- as_tibble(vignobles) %>%
  group_by(LIBCOM, NIVEAU) %>%
  summarise(N_Parcelles = n()
  ) %>%  
  pivot_wider(names_from = NIVEAU, values_from = N_Parcelles) %>%
  replace_na(list(`Bourgogne` = 0,
                  `Grand cru`=0,
                  `Premier cru`=0,
                  `Village`=0,
                  `Coteaux b.`=0))

## `summarise()` has grouped output by 'LIBCOM'. You can override using the
## `.groups` argument.

Jointure et proportions d’AOC

communes_final <- left_join(Vignobles_par_communes,
Vignobles_par_communes_par_AOC, by="LIBCOM") %>%
  left_join(communes,., by=c("nom_comm"="LIBCOM")) %>%
  mutate(share_Bourgogne = `Bourgogne` / N_Parcelles * 100,
         share_GdCru = `Grand cru` / N_Parcelles * 100,
         share_PmCru = `Premier cru` / N_Parcelles * 100,
         share_Village = `Village` / N_Parcelles * 100,
         share_Coteau = `Coteaux b.` / N_Parcelles * 100,
         share_PmGrCru = (`Premier cru`  + `Grand cru`) / N_Parcelles * 100,
         share_PmGrCruVill = (`Premier cru`  + `Grand cru` + `Village`) / N_Parcelles * 100)


head(communes_final)

# Prix moyen des vins estimés par commune
ggplot(data = communes_final) +
  geom_sf(aes(fill = Prix_moyen), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="red") +
  geom_text_repel(data = communes_final[communes_final$Prix_moyen > 100 ,],
                  aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Prix moyen par commune")

# Qualité moyenne des parcelles par commune
ggplot(data = communes_final) +
  geom_sf(aes(fill = MoyPond_b), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="goldenrod4") + 
  geom_text_repel(data = communes_final[communes_final$MoyPond_b > 80
                                      | communes_final$MoyPond_b < 65 ,], 
                  aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Qualités moyennes des vignobles par commune")

# Proportion de parcelles en grand cru par commune
ggplot(data = communes_final) +
  geom_sf(aes(fill = share_GdCru), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="goldenrod") +
  geom_text_repel(data = communes_final[communes_final$share_GdCru > 8 ,],
                  aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Part de grands crus par commune")

# Proportion de parcelles en premier ou grand cru par commune
ggplot(data = communes_final) +
  geom_sf(aes(fill = share_PmGrCru), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="gold") +
  geom_text_repel(data = communes_final[communes_final$share_PmGrCru > 30 ,],
                  aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")+
  labs(title = "Part de premiers et grands crus par commune")

Rapatriation des données communes au niveau vignoble

communes_final_tib <- as_tibble(communes_final) %>%
  select(-geometry)

vignobles_final <- left_join(vignobles,communes_final_tib, 
                              by=c("LIBCOM"= "nom_comm"))

summary(vignobles_final)

##     Concat             LIBCOM              NOM               NIVEAU         
##  Length:2391        Length:2391        Length:2391        Length:2391       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##     SURFACE          SCORE_b         Source             PrixMoy       
##  Min.   :0.0000   Min.   :18.40   Length:2391        Min.   :    5.0  
##  1st Qu.:0.1300   1st Qu.:67.90   Class :character   1st Qu.:   23.0  
##  Median :0.3100   Median :73.20   Mode  :character   Median :   36.0  
##  Mean   :0.4773   Mean   :71.89                      Mean   :  154.4  
##  3rd Qu.:0.6100   3rd Qu.:76.60                      3rd Qu.:   87.0  
##  Max.   :7.2000   Max.   :94.22                      Max.   :25960.0  
##                                                      NA's   :4        
##       gid          id_geofla       code_comm      insee_com   
##  Min.   :  563   Min.   :  563   412    : 158   21412  : 158  
##  1st Qu.: 6720   1st Qu.: 6720   295    : 128   21295  : 128  
##  Median :14392   Median :14392   054    : 126   21054  : 126  
##  Mean   :10484   Mean   :10484   464    : 121   21464  : 121  
##  3rd Qu.:14597   3rd Qu.:14597   541    : 121   21541  : 121  
##  Max.   :14749   Max.   :14749   (Other):1727   (Other):1727  
##  NA's   :10      NA's   :10      NA's   :  10   NA's   :  10  
##               statut       x_chf_lieu     y_chf_lieu      x_centroid  
##  Chef-lieu canton: 257   Min.   :8291   Min.   :66472   Min.   :8288  
##  Commune simple  :1998   1st Qu.:8342   1st Qu.:66555   1st Qu.:8338  
##  Sous-préfecture : 126   Median :8396   Median :66641   Median :8395  
##  NA's            :  10   Mean   :8405   Mean   :66654   Mean   :8404  
##                          3rd Qu.:8480   3rd Qu.:66758   3rd Qu.:8477  
##                          Max.   :8516   Max.   :66897   Max.   :8514  
##                          NA's   :10     NA's   :10      NA's   :10    
##    y_centroid       z_moyen        superficie     population       code_cant  
##  Min.   :66477   Min.   :222.0   Min.   :  90   Min.   : 0.200   05     :688  
##  1st Qu.:66553   1st Qu.:256.0   1st Qu.: 765   1st Qu.: 0.400   15     :515  
##  Median :66642   Median :300.0   Median :1006   Median : 0.700   23     :457  
##  Mean   :66655   Mean   :306.2   Mean   :1353   Mean   : 2.385   24     :374  
##  3rd Qu.:66751   3rd Qu.:355.0   3rd Qu.:1962   3rd Qu.: 1.600   06     :150  
##  Max.   :66899   Max.   :464.0   Max.   :3619   Max.   :21.900   (Other):197  
##  NA's   :10      NA's   :10      NA's   :10     NA's   :10       NA's   : 10  
##  code_arr    code_dept        nom_dept    code_reg        nom_region  
##  1   :1795   21  :2381   COTE-D'OR:2381   26  :2381   BOURGOGNE:2381  
##  2   : 586   NA's:  10   NA's     :  10   NA's:  10   NA's     :  10  
##  NA's:  10                                                            
##                                                                       
##                                                                       
##                                                                       
##                                                                       
##             Cote      aire_vignoble      MoyPond_b      N_Parcelles    
##  Côte D'or    : 152   Min.   :  5.25   Min.   :58.19   Min.   :  8.00  
##  Côte de Beane:1421   1st Qu.: 23.39   1st Qu.:69.98   1st Qu.: 63.00  
##  Côte de Nuit : 808   Median : 44.57   Median :71.39   Median : 88.00  
##  NA's         :  10   Mean   : 45.51   Mean   :71.07   Mean   : 93.89  
##                       3rd Qu.: 58.18   3rd Qu.:73.59   3rd Qu.:121.00  
##                       Max.   :101.94   Max.   :87.55   Max.   :158.00  
##                       NA's   :10       NA's   :10      NA's   :10      
##    Prix_moyen       Prix_median       Bourgogne       Grand cru    
##  Min.   :  16.55   Min.   :  17.0   Min.   : 2.00   Min.   :0.000  
##  1st Qu.:  26.11   1st Qu.:  26.0   1st Qu.:12.00   1st Qu.:0.000  
##  Median :  35.48   Median :  34.5   Median :21.00   Median :0.000  
##  Mean   : 154.16   Mean   : 112.9   Mean   :24.49   Mean   :1.351  
##  3rd Qu.: 140.79   3rd Qu.:  84.0   3rd Qu.:34.00   3rd Qu.:2.000  
##  Max.   :2151.34   Max.   :1447.0   Max.   :55.00   Max.   :9.000  
##  NA's   :10        NA's   :10       NA's   :10      NA's   :10     
##   Premier cru       Village        Coteaux b.    share_Bourgogne
##  Min.   : 0.00   Min.   : 2.00   Min.   : 0.00   Min.   :11.11  
##  1st Qu.: 8.00   1st Qu.:27.00   1st Qu.: 5.00   1st Qu.:19.05  
##  Median :17.00   Median :39.00   Median :10.00   Median :26.98  
##  Mean   :16.81   Mean   :38.15   Mean   :13.09   Mean   :26.08  
##  3rd Qu.:26.00   3rd Qu.:48.00   3rd Qu.:16.00   3rd Qu.:31.40  
##  Max.   :42.00   Max.   :62.00   Max.   :41.00   Max.   :75.00  
##  NA's   :10      NA's   :10      NA's   :10      NA's   :10     
##   share_GdCru      share_PmCru     share_Village    share_Coteau   
##  Min.   : 0.000   Min.   : 0.000   Min.   :11.54   Min.   : 0.000  
##  1st Qu.: 0.000   1st Qu.: 7.207   1st Qu.:32.77   1st Qu.: 4.688  
##  Median : 0.000   Median :19.672   Median :39.24   Median :11.570  
##  Mean   : 1.554   Mean   :17.136   Mean   :41.83   Mean   :13.398  
##  3rd Qu.: 2.500   3rd Qu.:24.324   3rd Qu.:48.98   3rd Qu.:17.857  
##  Max.   :11.111   Max.   :44.444   Max.   :76.54   Max.   :44.444  
##  NA's   :10       NA's   :10       NA's   :10      NA's   :10      
##  share_PmGrCru    share_PmGrCruVill          geometry   
##  Min.   : 0.000   Min.   :25.00     MULTIPOLYGON :2391  
##  1st Qu.: 7.207   1st Qu.:46.84     epsg:2154    :   0  
##  Median :23.529   Median :63.22     +proj=lcc ...:   0  
##  Mean   :18.690   Mean   :60.52                         
##  3rd Qu.:27.473   3rd Qu.:73.75                         
##  Max.   :55.556   Max.   :86.49                         
##  NA's   :10       NA's   :10

Modélisation multiniveaux

variable a expliquer: le prix moyen du 🍷

summary(vignobles_final$PrixMoy)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     5.0    23.0    36.0   154.4    87.0 25960.0       4

ggplot(data=vignobles_final, aes(x=PrixMoy)) +
  geom_histogram(fill="coral1")+
  labs(title = "Distribution des prix du vin")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 4 rows containing non-finite values (`stat_bin()`).

# Passage en log pour une distribution plus normale
vignobles_final$LogPrix <- log(vignobles_final$PrixMoy)

ggplot(data=vignobles_final, aes(x=LogPrix)) +
  geom_histogram(fill="coral3") +
  geom_vline(aes(xintercept=mean(LogPrix)),   
             color="black", linetype="dashed", size=1)+
  labs(title = "Distribution des log de prix moyens")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 4 rows containing non-finite values (`stat_bin()`).

## Warning: Removed 2391 rows containing missing values (`geom_vline()`).

variables explicatives niveau 1:

#Score qualité geographique (pente, geologie, pedologie, precipitations):
ggplot(data=vignobles_final, aes(x=SCORE_b)) +
  geom_histogram(fill="brown1")+
  labs(title = "Distribution des niveaux de qualité")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Surface:
ggplot(data=vignobles_final, aes(x=SURFACE)) +
  geom_histogram(fill="brown2")+
  labs(title = "Distribution des surfaces de vignoble")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Niveau d'aoc:
positions <- c("Coteaux b.", "Bourgogne", "Village",
               "Premier cru", "Grand cru")
ggplot(data=vignobles_final, aes(x=NIVEAU)) +
  geom_bar(fill="brown3") +
  scale_x_discrete(limits = positions)+
  labs(title = "Distribution des niveaux d'AOC")

#Source info prix:
ggplot(data=vignobles_final, aes(x=Source)) +
   geom_bar(fill="brown4") +
  labs(title = "Distribution des sources de données par vignoble")

variables explicatives niveau 2:

#Population de la commune
ggplot(data=communes_final, aes(x=population)) +
  geom_histogram(fill="aquamarine1", bins = 5)+
  labs(title = "Distribution des populations par commune")

#Part de grand et premier cru de la commune
ggplot(data=communes_final, aes(x=share_PmGrCru)) +
  geom_histogram(fill="aquamarine3", bins = 5)+
  labs(title = "Distribution des AOC premier et grand crus par commune")

#Cote de nuit ou cote de beaune:
pos_cote <- c("Côte de Nuit", "Côte de Beane", "Côte D'or")
ggplot(data=communes_final, aes(x=Cote)) +
  geom_bar(fill="aquamarine4") +
  scale_x_discrete(limits = pos_cote)+
  labs(title = "Distribution des communes par côte")

#Qualite moyenne de la commune
ggplot(data=communes_final, aes(x=MoyPond_b)) +
  geom_histogram(fill="darkslategrey", bins = 5) +
  labs(title = "Distribution des qualités moyennes")

Centrage-réduction-sélection des variables

vignobles_final_cr <- vignobles_final %>%
  mutate(LogPrix = scale(as.numeric(LogPrix)), 
         Surface = scale(as.numeric(SURFACE)), 
         Quality = scale(as.numeric(SCORE_b)), 
         AOC = fct_reorder(NIVEAU, SCORE_b),
         Cote = fct_reorder(Cote, SCORE_b),
         MeanQuality = scale(as.numeric(MoyPond_b)), 
         ShareGdCru = scale(as.numeric(share_GdCru))) %>%
  select(LogPrix, Surface, Quality, AOC,
         MeanQuality, ShareGdCru,
         LIBCOM, Source, Cote, Concat) %>%
  na.omit


summary_stat <- data.frame()

Modèle à un seul niveau

L1 <- lm(LogPrix~Quality+Surface+AOC+Source, data=vignobles_final_cr)

summary(L1)

## 
## Call:
## lm(formula = LogPrix ~ Quality + Surface + AOC + Source, data = vignobles_final_cr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0696 -0.6275 -0.2024  0.4016  3.9311 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.71893    0.09121  -7.882 4.88e-15 ***
## Quality             -0.01648    0.03248  -0.507 0.611999    
## Surface             -0.07580    0.01951  -3.885 0.000105 ***
## AOCBourgogne        -0.15699    0.07522  -2.087 0.036995 *  
## AOCVillage          -0.20275    0.08539  -2.374 0.017655 *  
## AOCPremier cru      -0.05133    0.11279  -0.455 0.649102    
## AOCGrand cru         2.33775    0.46750   5.001 6.14e-07 ***
## Sourceinterpolation  0.98998    0.05650  17.520  < 2e-16 ***
## Sourcewine_search    0.96382    0.46841   2.058 0.039735 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8659 on 2368 degrees of freedom
## Multiple R-squared:  0.2521, Adjusted R-squared:  0.2496 
## F-statistic:  99.8 on 8 and 2368 DF,  p-value: < 2.2e-16

Comme prévu, les vignobles de grandes surfaces ont des prix en moyenne moins élevés (“ce qui est rare est cher”), et les AOC “Grands Crus” sont significativement plus chers que les AOC Côteaux Bourguignons.
En revanche, la qualité physique des vignobles est non-significative dans le modèle des prix à un niveau, et les AOC “Bourgogne” et “Village” sont significativement moins chers que les AOC Côteaux Bourguignons.

Modèle nul à 2 niveaux

mnnull <- lmer(LogPrix ~ 1 + (1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mnnull)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: LogPrix ~ 1 + (1 | LIBCOM)
##    Data: vignobles_final_cr
## 
##      AIC      BIC   logLik deviance df.resid 
##   4487.1   4504.4  -2240.5   4481.1     2374 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4605 -0.3560  0.0610  0.4361  4.9734 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  LIBCOM   (Intercept) 0.8638   0.9294  
##  Residual             0.3611   0.6009  
## Number of obs: 2377, groups:  LIBCOM, 31
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)   0.1520     0.1677 30.9052   0.907    0.372

print(ranova(mnnull))

## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## LogPrix ~ (1 | LIBCOM)
##              npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>          3 -2240.5 4487.1                         
## (1 | LIBCOM)    2 -3371.2 6746.5 2261.4  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vrn2 <- summary(mnnull)$varcor[1][[1]][[1]]
vrn1 <- as.data.frame(summary(mnnull)$varcor)[2,]$vcov 
icc <- vrn2 / (vrn1 + vrn2)
icc

## [1] 0.7051989

La variance résiduelle du niveau commune est de vrn2 = 0.864.
La variance résiduelle du niveau vignoble est de vrn1 = 0.361.
Donc la correlation intraclasse ICC est de vrn2 / (vrn1 + vrn2) = 0.705.

Ou plus simplement:

performance::icc(mnnull)[[1]]

## [1] 0.7051989

Dans ce modèle vide à 2 niveaux, 29.5% de la variance est prise en charge par l’appartenance au niveau 2.

# Pour résumer:
summary_stat["MNNull","N"] <- dim(mnnull@frame)[1] 
summary_stat["MNNull","LogLik" ] <- summary(mnnull)$AIC[3][[1]] 
summary_stat["MNNull", "AIC" ] <-summary(mnnull)$AIC[1][[1]]
summary_stat["MNNull", "BIC" ] <- summary(mnnull)$AIC[2][[1]]
summary_stat["MNNull", "Deviance" ] <- summary(mnnull)$AIC[4][[1]]
summary_stat["MNNull", "Total_Var" ] <- vrn2 + vrn1
summary_stat["MNNull", "ICC" ] <- performance::icc(mnnull)[[1]]
summary_stat["MNNull", "InterClass" ] <- vrn1 / (vrn2 + vrn1)


summary_stat

Analyser les résidus de niveau 2

Extraction des résidus avec ranef, puis jointure et cartographie.

ranef <- ranef(mnnull)$LIBCOM %>%
  mutate(LIBCOM = as.character(rownames(.)),
         residuCom = as.numeric(`(Intercept)`))

mnnull_ranef <- left_join(communes_final, ranef, by=c("nom_comm"="LIBCOM"))

ggplot(data = mnnull_ranef) +
  geom_sf(aes(fill = residuCom), lwd=0.01, colour="white")  +
  scale_fill_gradient2(low = "#00008B", mid = "white", high = "#800020") +
  labs(title = "Résidus de niveau 2") +
  geom_text_repel(data = mnnull_ranef[abs(mnnull_ranef$residuCom) > 1,],
                        aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")

À comparer avec la carte originale des prix moyens par commune = différence de centrage-réduction.

ggplot(data = communes_final) +
  geom_sf(aes(fill = Prix_moyen), lwd=0.01, colour="white")  +
  scale_fill_gradient(low = "white", high="red") +
  labs(title = "Prix moyen du vin par commune") +
  geom_text_repel(data = communes_final[communes_final$Prix_moyen > 500 ,],
                        aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")

Modèle multiniveau avec variables de niveau 1 uniquement

mn1 <- lmer(LogPrix ~ Quality + Surface + Source + AOC +
              (1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mn1)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: LogPrix ~ Quality + Surface + Source + AOC + (1 | LIBCOM)
##    Data: vignobles_final_cr
## 
##      AIC      BIC   logLik deviance df.resid 
##   3038.4   3101.9  -1508.2   3016.4     2366 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4448 -0.5237 -0.0687  0.4612  4.6352 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  LIBCOM   (Intercept) 0.7839   0.8854  
##  Residual             0.1937   0.4401  
## Number of obs: 2377, groups:  LIBCOM, 31
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -0.50786    0.16599   36.31713  -3.060 0.004149 ** 
## Quality               -0.04720    0.01704 2347.41662  -2.770 0.005644 ** 
## Surface                0.04862    0.01028 2347.53697   4.728 2.40e-06 ***
## Sourceinterpolation    0.89024    0.02919 2346.53970  30.498  < 2e-16 ***
## Sourcewine_search      1.90564    0.24399 2347.88348   7.810 8.53e-15 ***
## AOCBourgogne          -0.07872    0.03881 2346.73517  -2.028 0.042646 *  
## AOCVillage            -0.17123    0.04409 2346.68467  -3.884 0.000106 ***
## AOCPremier cru        -0.07690    0.05793 2346.38460  -1.327 0.184518    
## AOCGrand cru           0.46195    0.24339 2347.85577   1.898 0.057822 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Qualty Surfac Srcntr Srcwn_ AOCBrg AOCVll AOCPrc
## Quality      0.152                                                 
## Surface     -0.010  0.229                                          
## Sorcntrpltn -0.168  0.072  0.212                                   
## Sorcwn_srch -0.006  0.060  0.191  0.094                            
## AOCBourgogn -0.197 -0.584 -0.073 -0.028 -0.028                     
## AOCVillage  -0.229 -0.720 -0.086  0.087 -0.022  0.817              
## AOCPremircr -0.237 -0.745 -0.096  0.234 -0.014  0.743  0.848       
## AOCGrandcru -0.067 -0.292 -0.258  0.002 -0.908  0.232  0.268  0.271

La relation entre prix et surface s’est inversée: les vignobles de grandes surfaces ont des prix en moyenne plus élevés.
Contre tout intuition, la qualité physique des vignobles est négativement (et significativement) associée aux prix dans ce modèle.
L’ajout des variables de niveau 1 a permis de réduire la variance totale de 20.2%.
La part de variance liée à l’appartenance communale est maintenant de 19.8%.

summary_stat["MNVar1Niv","N"] <- dim(mn1@frame)[1] 
summary_stat["MNVar1Niv","LogLik" ] <- summary(mn1)$AIC[3][[1]] 
summary_stat["MNVar1Niv", "AIC" ] <-summary(mn1)$AIC[1][[1]]
summary_stat["MNVar1Niv", "BIC" ] <- summary(mn1)$AIC[2][[1]]
summary_stat["MNVar1Niv", "Deviance" ] <- summary(mn1)$AIC[4][[1]]
summary_stat["MNVar1Niv", "Total_Var" ] <-
  summary(mn1)$varcor[1][[1]][[1]] + as.data.frame(summary(mn1)$varcor)[2,]$vcov 
summary_stat["MNVar1Niv", "ICC" ] <- performance::icc(mn1)[[1]]
summary_stat["MNVar1Niv", "InterClass" ] <- 1 -  performance::icc(mn1)[[1]]
summary_stat["MNVar1Niv", "PseudoR2" ] <- 
  (1 - (summary_stat["MNVar1Niv", "Total_Var" ]  / summary_stat["MNNull", "Total_Var" ] )) 
summary_stat

Résidus de niveau 1

mn1_residuals <- cbind(vignobles_final_cr, residuals(mn1))
colnames(mn1_residuals)[length(colnames(vignobles_final_cr))] <- "residusN1"

summary(mn1_residuals$residusN1)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.39613 -0.23049 -0.03025  0.00000  0.20294  2.03984

ggplot(data = mn1_residuals) +
  geom_sf(aes(fill = residusN1), lwd=0.01, colour="white")  +
  scale_fill_gradient2(low = "darkgreen", mid = "white", high = "darkorange") +
  labs(title = "Résidus de niveau 1") + 
  theme_minimal()

Résidus de niveau 2

ranef <- ranef(mn1)$LIBCOM %>%
  mutate(LIBCOM = as.character(rownames(.)),
         residuCom = as.numeric(`(Intercept)`))

mn1_ranef <- left_join(communes_final, ranef, by=c("nom_comm"="LIBCOM"))

ggplot(data = mn1_ranef) +
  geom_sf(aes(fill = residuCom), lwd=0.01, colour="white")  +
  scale_fill_gradient2(low = "#00008B", mid = "white", high = "#800020") +
  labs(title = "Résidus de niveau 2") +
  geom_text_repel(data = mn1_ranef[abs(mn1_ranef$residuCom) > 1,],
                        aes(label = nom_comm, geometry = geometry),
                  stat = "sf_coordinates",size=3, colour="black")

Modèle multiniveau avec variables de niveau 1 et 2

mn2 <- lmer(LogPrix ~ Quality + Surface + Source + AOC +
               MeanQuality + ShareGdCru + Cote +
              (1 | LIBCOM), data=vignobles_final_cr, REML=F, na.action=na.omit)
summary(mn2)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## LogPrix ~ Quality + Surface + Source + AOC + MeanQuality + ShareGdCru +  
##     Cote + (1 | LIBCOM)
##    Data: vignobles_final_cr
## 
##      AIC      BIC   logLik deviance df.resid 
##   2999.6   3086.2  -1484.8   2969.6     2362 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4420 -0.5218 -0.0649  0.4570  4.6376 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  LIBCOM   (Intercept) 0.1684   0.4104  
##  Residual             0.1937   0.4401  
## Number of obs: 2377, groups:  LIBCOM, 31
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -1.05896    0.11393   42.28250  -9.295 8.91e-12 ***
## Quality               -0.04788    0.01704 2346.30069  -2.809 0.005004 ** 
## Surface                0.04860    0.01028 2350.17563   4.727 2.41e-06 ***
## Sourceinterpolation    0.89129    0.02919 2347.54433  30.537  < 2e-16 ***
## Sourcewine_search      1.88415    0.24384 2353.35431   7.727 1.62e-14 ***
## AOCBourgogne          -0.07774    0.03880 2347.94052  -2.003 0.045241 *  
## AOCVillage            -0.17095    0.04408 2347.70878  -3.878 0.000108 ***
## AOCPremier cru        -0.07587    0.05793 2347.03637  -1.310 0.190419    
## AOCGrand cru           0.47943    0.24325 2352.82532   1.971 0.048845 *  
## MeanQuality           -0.04275    0.08039   30.57464  -0.532 0.598741    
## ShareGdCru             0.36251    0.08666   30.27123   4.183 0.000227 ***
## CoteCôte D'or          0.62786    0.26975   31.72140   2.328 0.026477 *  
## CoteCôte de Nuit       1.14890    0.17176   30.09307   6.689 2.04e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

La relation entre prix et surface est toujours inversée: les vignobles de grandes surfaces ont des prix en moyenne significativement plus élevés.
La qualité physique des vignobles est toujours négativement (et significativement) associée aux prix dans ce modèle.
L’AOC “Grand Cru” est positivement (et significativement) associée aux prix du vin dans ce modèle, mais moins fortement que la part de grands crus présents dans la commune (ShareGdCru).
Par rapport à la Côte de Beaune, les vins de la Côte de Nuit se vendent très significativement plus chers.
L’ajout des variables de niveaux 1 et 2 a permis de réduire la variance totale de 70.4%.
La variance interclasse est maintenant de 53.5%.

summary_stat["MNVar2Niv","N"] <- dim(mn2@frame)[1] 
summary_stat["MNVar2Niv","LogLik" ] <- summary(mn2)$AIC[3][[1]] 
summary_stat["MNVar2Niv", "AIC" ] <-summary(mn2)$AIC[1][[1]]
summary_stat["MNVar2Niv", "BIC" ] <- summary(mn2)$AIC[2][[1]]
summary_stat["MNVar2Niv", "Deviance" ] <- summary(mn2)$AIC[4][[1]]
summary_stat["MNVar2Niv", "Total_Var" ] <-
  summary(mn2)$varcor[1][[1]][[1]] + as.data.frame(summary(mn2)$varcor)[2,]$vcov 
summary_stat["MNVar2Niv", "ICC" ] <- performance::icc(mn2)[[1]]
summary_stat["MNVar2Niv", "InterClass" ] <- 1 -  performance::icc(mn2)[[1]]
summary_stat["MNVar2Niv", "PseudoR2" ] <- 
  (1 - (summary_stat["MNVar2Niv", "Total_Var" ]  / summary_stat["MNNull", "Total_Var" ] )) 
summary_stat